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Abstract 

Hot spots in tumors are regions of high vascular density in the center of the tumor and their 
analysis is an important diagnostic tool in cancer treatment. We present a model for vascular 
remodeling in tumors predicting that the formation of hot spots correlates with local 
inhomogeneities of the original arterio-venous vasculature of the healthy tissue. Probable 
locations for hot spots in the late stages of the tumor are locations of increased blood pressure 
gradients. The developing tumor vasculature is non- hierarchical but still complex displaying 
algebraically decaying density distributions. 
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1 Introduction 



Tumors that grow inside vascularized tissues remodel their vascular environment in a characteristic 

way. Tumor cells secret growth factors, the most prominent being vascular endothelial growth 
factor (VEGF) that diffuse into the surrounding tissue. This initiates a vascularization program 
that is denoted as tumor- induced angiogenesis (Carmeliet and Jain, 2000). The lumen of blood 
vessels is surrounded by endothelial cells (ECs), which express VEGF receptors that stimulate 
EC migration and proliferation upon binding to VEGF. Existing blood vessels generate sprouts 
that migrate through the extracellular matrix, along a chemotactic gradient if present or guided 
by receptors like those from the Eph family and their ephrin ligands (KuUandcr and Klein, 2002; 
Carmeliet and Tessier-Lavigne, 2005) or just randomly in the absence of any guidance clue. As 
a result new blood vessels are formed in the peritumoral region, where the microvascular density 
(MVD) is increased. This improves the oxygen and nutrient supply of the periphery and leads to 
proliferation and expansion of the tumor. 

Inside the tumor the formation of new vessels deceases in spite of an abundance of growth 
factors and the vascularization program is switched from sprouting angiogenesis to circumferential 
vessel growth. Experiments demonstrate that certain guidance molecules like EphB4 expressed by 
tumor ECs act as a negative regulator of blood vessel branching and vascular network formation 
(Erber et al., 2006). In addition to a drastic dilation of the blood vessels in the tumor center 
a massive vessel regression inside the tumor is observed (Holash et al, 1999; Maisonpierre et al., 
1999) leading to a much lower MVD in the tumor center. Regression of blood vessels inside the 
tumor can be evoked by the presence of angiogenic agonists like angiopoetin (Holash 1999a, 1999b) 
or by mechanical solid stress in conjunction with a leaky and disorganized vessel structure leading 
to vessel collapse. 

This scenario deviates from the prevailing picture that most tumors and metastases originate 
as small avascular masses that belatedly induce the development of new blood vessels once they 
grow to a few millimeters in size (Folkman, 1971; 1990). Instead tumors rapidly coopt existing 
host vessels to form an initially well-vascularized tumor mass, which in later stages attains its 
characteristic morphology, which is characterized by (Paku, 1998; Holash 1999a, 1999b; Dome et 
al., 2002): The tumor is compartmentalized into a perimeter region with MVD much higher than 
in the normal tissue, a well vascularized tumor periphery, and a tumor center with a very low 
MVD. In the tumor center necrotic regions are threaded by a few very thick vessels surrounded 
by 100-200 /Ltm thick cuffs of viable tumor cells. This is expected to be the general course of 
development of a tumor vessel network that grow in a vascularized tissue. 

Recently two of us introduced a mathematical model (Bartha and Rieger, 2006) predicting the 
dynamics of vascular network formation based on the mechanisms identified experimentally and 
described in part above. Its basic ingredients are a dynamic network, in which links, representing 
vessels, can be added or removed or altered (in diameter and other parameters), a cellular au- 
tomaton representing the growing tumor, and two concentration fields mediating the interaction 
between the network and the tumor cells. The links of the network, modeled as pipes carrying 
an ideal pipe flow, constitute the sources of the oxygen or nutrient concentration field, and the 
tumor cells are sources of the growth factor concentration field. Generation and removal of vessels 
or tumor cells are stochastic processes with transition probabilities that depend on local concen- 
trations of oxygen and growth factor, blood flow and duration of under-oxygenation. The model 
displays the characteristic features of tumor vasculature described above and it is robust against 
parameter variations and various refinements (Welter et al., 2007) as well as the choice of space 
dimensionality (Lee et al., 2006). 

Here, the global characteristics of the emerging morphology of the tumor vasculature, like the 
radial dependencies of MVD, tumor density, blood flow etc., were independent of the details of 
the original vasculature. However, local characteristics, like the formation of highly vascularized 
regions in the tumor center, the so called "hot spots" (Weidncr, 1995; Bclien et al., 1999), could not 
be analyzed since the original vessel network was assumed to form a regular lattice and the original 
blood flow imprinted via boundary conditions was directed along the diagonal of the underlying 
lattice. As a result, the surviving thick vessels threading the fully developed tumor where also 
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running along this diagonal and these represented the only hot spot, always occurring in the same 
location of the tumor network. Since the analysis of hot spots are an important diagnostic tool in 
cancer treatment (Weidner, 1995; Belien at al., 1999; Pahernik et al., 2001) it is highly desirable 
to understand the mechanisms determining their development in growing tumors. 

The main hypothesis we intend to push forward in this work is that hot spot formation in 
tumors is connected to local characteristics of the original vascular network. For this reason 
it appears mandatory to start with a realistic initial vasculature, reflecting the; typical arterio- 
venous organization of the blood vessel network in healthy tissues. The construction of such a 
network consisting of two interdigitating hierarchical trees, an arterial and a venous tree, under 
the constraint that the lowest level, the capillaries, provide a homogeneous supply of oxygen in 
the tissue, is a highly non-trivial task already in a two-dimensional set-up. A mathematical model 
for the stochastic generation of an arterio-venous network was developed in (Godde and Kurz, 
2001) and a modified version of it is used as the starting configuration of the vessel network in 
our model. This is based on the ideas developed in our previous work (Welter et al, 2007) for 
a regular starting network and adapted to the presence of arteries and veins with varying vessel 
radius and blood flow parameters. The purpose of this work is to describe the transformation of 
the initial arterio-venous vasculature into a non-hicrachical tumor vasculature with that displays 
still complex properties characterized by algebraically decaying density distributions. In particular 
it leads to the formation of hot spots and the most probable location of the regions with highest 
vascular density can be predicted from an analysis of pressure gradients in the original network. 

The paper is organized as follows: In the next section the model is defined, including the method 
to construct an initial arterio-venous network. Section 3 presents the results for global properties 
of the tumor vasculature for a choice of parameters that is guided by experimental data for 
melanoma. These results include a discussion of the emerging morphologies; radial distributions of 
vessel density, vessel radius, tumor density, flow rates, shear forces etc., vessel statistics, parameter 
dependencies and drug flow. In section 4 we analyze the spatial inhomogeneities and the hot spot 
formation of the tumor vasculature. Section 5 concludes the paper with a discussion. 

2 Mathematical Model 

Our mathematical model contains discrete elements and continuous fields. The tumor is repre- 
sented by a set of occupied sites on a lattice and the vascular system by a network of connected 
tubes which occupy lattice bonds. The continuous part involves a O2 field as nutrient supply and a 
growth factor (GF) field for angiogenic signaling. The temporal evolution is governed by stochastic 
processes according to which the automaton state is siiccessively updated. An exhaustive descrip- 
tion is given in our previous publication (Welter et al., 2007), but since the introduction of an 
arterio-venous network and some refinements have lead to a few changes we give a comprehensive 
overview here. 

2.1 Definition of the system state 

As illustrated in Fig.l, we study a quasi two-dimensional case, where the arrangement of vessels 
and tumor cells (TCs) is subject to a regular triangular lattice with lattice-constant 6 and size 
I X I = n. Each lattice site can be occupied each with a single tumor cell. We define T as the set 
of sites currently occupied by viable TCs. Non-occupied sites represent normal tissue, which we 
assume to be homogeneous although this is not the case in reality. Likewise lattice bonds can be 
occupied with vessels. The vasculature is modeled as network of connected ideal tubes. We can 
formally describe its topology by a graph G = {N,V), where the edges mV = {{i,j), i,j G N} 
represent vessel segments which we also denote as vessels. The elements in N refer to attached 
junction nodes. Note that vessels run parallel to the lattice axes and that node locations coincide 
with lattice sites. We further refer to various properties via sub-scripts, for example the site 
location of a node i is given by a;,, while qaifa^la denote flow-rate, radius and length of some vessel 
a. 
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Blood flow is an important part of our modeling. Wc assume steady-state laminar flow through 
the network. Given the pressure gradient Ap = {pi — pj)/l between the ends of a vessel 
Hagen-Poiseuille's Law determines the flow rate q oc r'^Ap/rj and wall shear-stress fa oc rAp. To 
account for the complex non-Newtonian fluid behavior of blood it is common to introduce an 
efi'ective viscosity rj = rjpiasma VreiiH, r) depending on local vessel radius r and hematocrit H. We 
use the empirical formula derived by Pries et al. (1994), where the precise deflnition can be found. 
For simplicity we neglect the phase separation efifect and set the hematocrit to 0.45, the average 
in the human body. To find the nodal pressures we solve the system of linear equations which 
arises from the application of mass-conservation at junctions. Appropriate boundary conditions 
must be present. Our arterio-venous trees require them only at the ends of the major supply-and 
drain vessels. We have chosen to set the respective nodal pressures to fixed values P^''\r) based 
on the vessel radius. 

Vasculature and tumor interact via the growth factor field Cg and the O2 field Cq. The latter is 
given by the reaction diffusion equation Aco — kCq + a{co — Co) = 0, where sources are present 
at vessel locations and sinks are distributed everywhere. We solve for steady states, because the 
relaxation time is of the order of seconds while configuration changes at the cellular level take 
hours. The coefficients have the following meaning: k is a tissue dependent consumption rate, 
Co the blood oxygen level and a the source strength coefficient incorporating local vascular 
area-fraction and wall permeability. Thereby a number of assumptions were made: (i) O2 uptake 
is supply-limited so that a linear approximation to the I0W-PO2 regime of a Michaelis-Mentcn 
relationship is reasonable, (ii) Blood-P02 is constant throughout the vasculature, (iii) Vessels 
contribute equally to O2 supply, i.e. a is a constant a^°' > at occupied sites, ignoring differences 
in wall thickness, surface area, etc. Note that the presence of the sink-term kCq (k > 0) would 
lead to exponentially decaying distributions around point sources. Indeed, experiments show an 
approximately linear PO2 decay from the outer vessel wall ranging ca. 150 jim into the tissue. 

In analogy, we assume that due to binding and degradation GF has a limited range which is 
explicitly given by i?^^^ . Let us denote T„o as the set of underoxygenized sites among x GT. These 
are characterized by Co{x) < 6^^°^^ where 9^™''^ is a threshold parameter. Under the assumption 
that TCs produce GF at a constant rate the formulation of a diffusion-reaction equation leads to 
a simple solution via a Greens function approach where Cg is given as Cg{x) = J2x'eT di""" ~ ^')- 
For simplicity, convenience and eflaciency, we use a linearly decaying "Greens function": g{r) oc 
max[0, l-r/ii(s)]. 

2.2 Definition of the dynamical processes 

Tumor and vascular network dynamics are governed by the stochastic processes described below. 
In practice, when we analyze the model by Monte-Carlo simulations, we iteratively update the 
system by a fixed time step At = Ih. Per step, a sweep is done for each process. After that the 
configuration has changed and therefore O2 and GF fields are recomputed. This is reasonable if 5 
is sufficiently smaller than the rates at which changes are introduced. 

TC Proliferation: New TCs occupy empty neighbor sites x with probability Af/i^™'^ if 

the local oxygen level is high enough Co{x) > 6^^°^\ This resembles the behavior of real 

tumors (Bru et al., 2003; Drasdo and Hohmc, 2005) where proliferation is restricted to a 
small band behind the invasive edge due to space and nutrient constraints. {Fig. 2a) 

TC Death: TCs are removed with probability p^'^"*'*'' ~ 1/2 if the local O 2 level is less than 
0(deaih) longer than t^r^Q ■ Since TCs adapt to low oxygen conditions (Iyer et al., 1998), 
is very small. For simplicity the survival time t'^j^Q has been given a 
fixed value, in spite of strong variations among genotypes (Yu et al., 2002). {Fig. 2b) 

Angiogenesis: The presence of GF stimulates angiogenic sprouting at a vessel occupied site x 
in the following way. First it is required that the local GF concentration Cg{x) is larger 
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than the threshold 9g° . Evidently no sprouting happens inside the tumor (Holash et al. , 
1999). Indeed molecular pathways (via EphB4 and its ligand ephrinB2, expressed in tumor 
endothelial cells) have been identified recently which are related to the switch from sprouting 
to circumferential growth (Erber et al., 2006). Therefore it is also required that parent 
vessels have not been inside the tumor for longer than t^™*'^'*'' . We characterize a position 
as "inside" if radial distance to the tumor center is smaller than the extent of the tumor in 
the respective direction. A vessel is inside the tumor if the midpoints of at least 50% of all 
occupied bond are inside. Our vessels carry a type information: they are either artery, vein or 
capillary. Vessels created via sprouting as described here are always of capillary type, while 
the type of the initial vessels is given by construction of the arterio-venous trees. For our 
base case, we require the parent vessel to be either of venous or capillary type, automatically 
including active (capillary-)sprouts. These spouts can thus already branch out into sub- 
sprouts. We exclude arteries because the thicker endothelium and presence of smooth-muscle 
cells makes sprouting less likely. The last requirement is that adjacent bifurcations must be 
at least /(^p*") ^vh away since endothelial cells close to bifurcations also do not form sprouts. 
If those criteria are fulfilled an initial segment of length 5 and radius r^""*) is appended with 
probability At/t^^^°^*\ possibly splitting the parent vessel. At successive iterations more 
such segments are added to the tip with probability At/t^^^°^*K On contact, it connects to 
other vessels forming a potentially blood perfused loop. Thereby the initial sprout direction 
is given by the steepest GF descent (Gerhardt et al., 2003). Furthermore, sprouts cannot 
extend indefinitely (Nehls et al., 1998). Therefore a "countdown" of s^™'^^-' hours terminates 
the sprouting behavior upon expiration (see below). The respective internal timer variable 
is inherited from the previous sprouting tip. Note that tumors grow by co-opting vessels but 
no explicit modeling is required for this since TCs proliferate in a separate layer without 
mechanical interaction. (Fig. 2c, d) 

Vessel Collapse: Vessel in normal arterio-venous networks have differing amounts of support 
structures around the endothelium: the basement membrane, pericytes and smooth-muscle 
cells. In tumors interaction/adhesion between these components is often disrupted (McDon- 
ald and Choyke, 2003). There one finds vessels with fragmented and multi-layered basement 
membranes and detached pericytes. Further the lumen of such vessels is often completely 
collapsed and/or the endothelium is also found in a state of regression. Therefore, we believe 
that it is physiologically sound that the survival time of initially healthy vessels depends on 
how well the vessel wall is developed, i.e. capillaries with just the membrane around the 
endothelium should collapse earlier than thicker arterioles. We realize this by a "degree- 
of-maturation" variable w which is initialized with data for the wall thickness in normal 
vessels (Pries et al., 2005), depending on vessel radius r. Using a rough approximation, we 
set w = 2r(0.65 — 0.21og(2r)). Inside the tumor, vessels are continuously degraded at the 
rate Aw. While w is larger than the threshold it;^'^\ a vessel does not regress, meaning the 
segment is not removed from the network. Furthermore, physiological wall shear stress levels 
can dose dependently inhibit endothelial cell apoptosis, while long term reduction can cause 
vessel regression (Dimnieler and Zeiher, 2000). We include this by constraining removals to 
vessels with shear forces / less than the threshold /*^^^. In tumors, impaired blood flow and 
thereby low shear stresses might be the result of increased solid pressure compressing the 
vessels. Also our previous work has shown that shear force correlated collapses lead to real- 
istic network morphologies. If collapse inhibiting effects do not apply, segments are removed 
with probability . For a dose dependent effect we set p^'^-' = p^^^n = (1.0 — ///'^^■'), 

with the constant parameter ]\[Q^g that sprouts are completely excluded from the 

above processes. In section 13.31 we also discuss the effect of restricting collapses to a thin 
band behind the invasive edge, using the collapse probability as defined there. (Fig.2e) 

Vessel Dilation: Exposure to growth factors is also associated with unphysiological vessel dila- 
tion, in particular in tumors which is related to the mentioned switching to circumferential 
growth inside the tumor (Erber et al., 2006). Experiments (Dome et al., 2002) indicate that 
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vessel diameters increase continually to an upper limit. In the model the dilation process 
works as follows: Each bond occupied by a non sprouting vessel a increases the radius ra 
by Ava with probability At/t^^Q^ if ra < r^"^°'^\ if the center of the bond is inside the tu- 
mor and if the local GF concentration is larger than 6^^°''\ The increment Ara = jl-nla 
corresponds to the surface area of an additional EC contributed to the total surface area of 
the segment. To account for surface tension, a smoothing effect is generated by dilating the 
thinnest vessel at junction bonds. {Fig.2f) 

2.3 Arterio- venous tree generation 

Vasculature in living tissue exhibits a tree like structure. Few thick arteries branch out into 
arteriolar microvessels. Terminal branches are connected to the capillary bed, a dense network 
consisting of thin vessels where most of the exchange with the surrounding tissue happens. Further 
upstream blood is collected in venoulcs which fuse into thick veins. In the sense of some optimality 
measure, the design goal of such structure is to provide sufficient amount of nutrients to tissue 
jwhile minimizing the effort to keep the circulatory system operating. 

Numerous ways to generate vasculature in silico have been proposed in the litca'ature for exam- 
ple fractal approaches based on L-Systems (Mandelbrot, 1991) or optimization based geometrical 
construction (Schrener and Buxbaum, 1993). Most of which concentrate on building the supply- 
ing arteries only, neglecting explicit modeling of capillaries and veins. A deterministic manual 
construction of a single configuration would have been unsuitable and difficult with respect to the 
desire to obtain the flexibility to create a wide range of configurations while retaining physiological 
properties. 

Godde and Kurz, (2001) presented a method to grow and remodel vascular trees stochastically 
according to probability functions that depend on local system properties. The simultaneous 
construction of arterial and venous trees makes their approach well suited for generating an initial 
vasculature for our tumor simulator, where a complete network including draining venoules is 
required. 

The general idea is to grow purely random arterial- and veneous trees on a lattice, under 
exclusion of already occupied sites. Followed by iterative shear-stress guided further growth and 
regression at the tree leafs, whereby a well perfused space filling network evolves. 

2.3.1 Network construction 

The mathematical model for the vasculature is identical to the definition in the tumor simulation: 
the network is identified with a graph plus additional geometric and hemodynamic properties, 
embedded in a triangular lattice. 

Our initial condition for the growth stage consists of several prescribed segment chains for 
arteries and veins which become the major vessels in the upper tree levels. Their length and 
position is drawn from a uniform distribution within reasonable bounds. 

The basic structural element for further growth is a tripod consisting of three vessel segments 
arranged in 120° angles. These tripods are appended successively at randomly chosen tree leafs. 
Exceptions are the initial chains where all included sites are candidates for adding elements. 
Thereby overlap with already occupied sites/bonds is forbidden, and so the process continues 
until there are no more valid configurations where tripods could be added. The resulting network 
consists of at least two binary trees with at least one arterial and one venous side. 

To compute blood flow vessel radii must be known. 'Murray's Law' (Murray, 1926) relates the 
radius of a parent vessel a to the radi of branches b,c. a" = + c", where a = 2.7 is a realistic 
value as used by Godde and Kurz. The radii of the tree leafs arc preset to 4 /im for arteries, and 5 
fjim for veins. Thus all other radii can be computed in a recursive algorithm which visits a vessels 
once its descendants have been processed. 

To get a complete network, capillary segments have to be added. We simply loop through the 
bonds on the lattice and add a connection a with = 3.5 /um if it connects an arterial tree with 
a venous tree. As additional constraint a connection is not added if more than three vessels would 
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meet at a site. Capillaries are not constraint to connect leafs nodes exclusively, however all vessels 
at the respective node pair must fulfill r < r^^"^*). 

The second stage of construction consists of iterating: radii determination, capillary creation, 
flow rate and shear force computation, capillary removal, and remodeling sweep, until the system 
reaches a steady state. Our remodeling procedure works as follows: One of growth, death, idle 
is drawn for each segment according to the respective probabilities pg,pd,l — Pg — Pd, where 
Pg+Pd < 1- Relating these to shear stress / such that Pg y f and gd\ f leads to expansion in well 
perfused branches while non-perfused branches regress and make room for growth of interdigitating 
patterns. The probability relationship we used here is as follows: 

I ^ ln/-mina{ln fa} 

maXa {In /a } — miiia {In /a } 

Pg ^if + Po)-r 

9d ^(l-.f+Po) 

Pg ^ P~g/{Pg+P~d) 

Pd ^ Pg/{Pd+Pd) 



To prevent unnatural situations where thick vessels have too many capillary branches while still 
obtaining homogeneous MVD there is the radius term f which decays from 1 at r = 4 to at 10 
/xm. The offset po is added for two reasons: to lower the rate at which un-circulated initial vessels 
regress, and to add fluctuations to reduce deadlocks in high shear stress regions. Leaf segments 
labeled as dying are removed. The removal is constraint to leafs in order to maintain the binary 
tree structure. Next, new elements are added to growing segments, whereby one of all admissible 
configurations is randomly picked. Beside tripods we also allow single segments to be added in 
order to allow vascularization of regions where bottlenecks in between major vessels would block 
access. 

2.3.2 Construction results 

In order to study typical tumor-network morphologies we must simulate a sufficiently large domain: 
12x10 mm. To obtain a specific MVD we adjust the tree generators bond-length (jfs'^") and set 
the lattice size accordingly. Due to the two-dimensionality of our system we measure MVD as 
the fraction of occupied sites. Furthermore we chose 6^^^'^^ to be a multiple of the tumor lattice 
bond-length S so that the network can be taken to the tumor model by superimposing the lattices. 
Based on the average for human skin, 100 vessels per mm^ cross-section (Dome et al., 2002) which 
implies 10 /zm inter-vessel spacing, we set 5^^^"'^ = 60 /xm and I = 200. 

Due to the hierarchical design of real vasculature, blood-pressure correlates well with vessel 
radius. Plotted, it exhibits a sigmoidal shape (Pries et al., 1995). In order to prevent physiologi- 
cally inconsistent pressure gradients across our networks, we adapt the boundary pressure by an 
approximate formula. (p(r) = 0.133(18 -|- 72/ (1 -|- exp((r -|- 21/im)/16/im))) kPa; r is negated if it 
is the radius of an arterial segment). 

As illustrated in Fig. 3, we have chosen four basic layouts. The inlet schematics indicate bounds 
for the positions and lengths of the initial chains. These are drawn from uniform distributions 
within a prescribed range, whereas their types (arterial or venous) are unalterable. The initial 
conditions were deliberately chosen, motivated by the observation that vasculatures exhibit rel- 
atively thick straight vessels down to a certain level in the hierarchy as can be seen in photos 
from the CAM (Mironov et al., 1998), or in recent 3d vascular imaging experiments (Cassot et al., 
2006). 

Quantitatively we have compared our data with scatter plots of hemodynamic properties shown 
in (Godde and Kurz, 2001). They are near identical to the original and therefore not shown directly. 
Instead we refer to Fig. 10 where tumor vessels are also included. 
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2.4 Parameters 



In what follows we analyze the model with parameters guided by experimental data (Dome et al., 
2002) for human malignant melanomas. 

The lattice bond length S is set to 10 /im, the typical diameter of TCs and ECs. The initial 
tumor is an cdcn-grown cluster of 1000 colls and 300 /im in diameter. The initial vasculature has 
been discussed. Growth and death rate parameters as well as oxygen and growth factor diffusion 
parameters have been discussed in detail in (Welter et al., 2007). 

To recapitulate, we set the TC proliferation time to 10 h, sprout initiation/extension 

time f^p'^""*-' to 5 h, the time until regression of loose sprouts gC™"^) to 100 h (Nchls et al. 1998), 
the minimum distance between bifurcations i to 20 /xm, the initial sprout radius r*^™**) to 4 
/im, the time until sprouting switches to circumferential growth f^^**'^'*-' to 24 h, EC proliferation 

time for circumferential growth t^^Q^ to 40 h and the maximum radius r^™"^) to 25 fxm. For the 
continuous fields we set the O2 source coefficient a*-'^-' to 0.002, the consumption coefficient n to 
^{Nmm) _ ^gg ^-2 normal tissue and to k*^"^-* = 4k(^™™) in tumor tissue to get halve the 

diffusion distance. Here the blood oxygen level Co can be chosen arbitrarily because it appears 
as global scaling factor in the solution. The model however can be adjusted accordingly by scaling 
the threshold parameters. Therefore we set it as previously to the hematocrit value Co =0.45 
(which is constant here). We set the threshold for TC proliferation 6^^°''^ to 0.3 ~ 0.9(co) and 
threshold for extreme under-oxygenation ^('^^°*'') much smaller to 6^^°''^ /lO. The TC survival 
time under hypoxia is set to 100 h. The growth factor radius i?^^^ is set to 200 fj,m, and the 

vessel proliferation threshold Og^™''^ is set very low to 10~^ so that all vessels within the full GF 
radius are affected. Parameters related to vessel regression are as follows. The critical shear force 

average shear stress is (/) « 10 Pa. The collapse probability 
p{c,max) jg relatively high to 0.1. The dematuration rate Aw is set to 0.04 /im/h, resulting in 
regression delays from 93h (4/im vessels) to 625h (50/um vessels). 



3 Results 

We run 10 Monte-Carla simulations for each of the four initial vascular setups as explained above 
and shown in Fig. 3. We generated different initial networks for each run by changing the random 
number seed for the construction algorithm. The vessel density distribution in the vascularized 
regions is generally very homogeneous by design since we attempt to maximize the lattice occu- 
pation. Normally the networks fill less than the entire rectangular domains. We accept this if the 
support area is sufficient for tumor growth till the end at t =1200h 

Fig. 4 shows the evolution of the model using vasculature configuration (a). Not shown is 
the very first stage of tumor growth which includes the following. High tumor O2 consumption 
leads to decreased O2 levels inside the nucleus and consequently enables vascular remodeling via 
growth factor production. Angiogenic sprouting increases the MVD around the tumor, while after 
a few hours of prolonged exposure to GFs, central capillaries show noticeable increased diameters. 
The first image shows the tumor at t =200h where a dense capillary plexus has been created by 
angiogenic sprouting around the tumor nucleus. Vessel collapses are not yet occurring there. By 
t =100h the wall maturation w decreased sufficiently though to allows first collapses of the thinnest 
capillaries with inadequate shear force. We observe this in other simulation runs. The process of 
angiogenesis and collapse proceeds as the tumor expands. At t =400h first regions are visible where 
TCs have died due to hypoxia. Thick arterioles and venoules don't yet collapse due to u> > w^^^ 
though they may have low shear force levels. Fig. 5 shows the final tumor at t =1200h and Fig. 6 
shows tumors of the other network configurations. Fig. 7 shows a magnified view on the boundary 
region of the tumor in Fig. 5. The tumor masses expand approximately disc shaped. Because 
of the simple oxygen distribution model where all vessels release O2 homogeneously, this fulfills 
our expectations. Normally nutrient exchange happens at the thinnest capillaries. Venoules and 
arterioles would naturally not contribute to tissue supply as much. This however is not practical 



7 



in a two dimensional model because non-oxygen-releasing vessels would be impassable barriers to 
the tumor. Therefore we require that O2 content, and wall-permeability is equal for all vessel. The 
resulting network morphologies show the typical high MVD periphery /peritumoral region and low 
density center. A general feature of our model seems to be the formation of arterio-veneous shunts 
cause by dilation of vascular pathways that connect to initially available arterioles/ venoules. Those 
paths can consist of neo-vasculature as well as initial capillaries and other vessels. As we expect, 
most vessels form isolated straight threads, however regions where vessels are densely clustered 
are also observable. 

3.1 Radial distributions 

For comparison Fig. 9 shows quantities averaged within concentric shells emanating from the tu- 
mor center. MVD and TC density arc computed as the fraction of occupied sites within a shell. 
Hemodynamic properties are averaged over occupied bonds only. Mean values of the microvas- 
cular density MVD and vessel radius r, reminiscent of results from earlier work, agree well with 
experiments (Dome ct al., 2002) as we have discussed for similar data in (Bartha and Riegcr, 

2006) . Dome et al. measured the MVD and MVR in three distinct regions: the central region, 
a 100/xm wide peripheral band just behind the invasive edge, a 200/im wide peritumoral region 
outside the invasive edge. In the central region, they found 25% MVD of normal tissue, and up 
to 200% in the peritumoral region. The vessel perimeter grows linearly from 50/xto and assumes 
a plateau at 200/xm by day 15. 

In contrast to results for a regular vascular network (Bartha and Riegcr, 2006; Welter ct al., 

2007) , flow rates and shear force now show a plateau like the vessel radius which is expected since 
the flow boundary condition on the regular network lead to unrealistic star shaped morphologies, 
directing all blood flow through the center. The hierarchical networks here do not introduce 
such obvious problems. Not shown is data for the blood-pressure gradient dp/dl which decreases 
monotonically toward the center by more than one order of magnitude. The blood flow rate q is 
proportional to r'^dp/dl and increases on average toward the center beyond a small local minimum 
at the invasive edge. Thus the dependency of the radius outweighs the pressure drop. The 
shear force dependency is r dp/dl, leading to a decrease respectively. 

3.2 Vessel Statistics 

Fig. 10 shows scatter plots of hemodynamic variables against the vessel radius r. While providing 
more information about their distribution over vessels, such plots are also common in the literature. 
Compared to (Goddc and Kurz, 2001) we get a good match. The data is from a single simulation 
run. Data from the initial vasculature and from the tumor vasculature at t =1000h are both 
displayed in the same plots. 

For the initial network we generally observe that the variance of the flow related parameters 
increases drastically towards the capillaries. While this might be an artifact of the artificially 
constructed network, averaged values show physiologically sound characteristics. 

The dilation of tumor internal vessels leads to samples clustered at the maximum dilation radius 
^^{max) _ 25yrini) and therefore one observes larger ranges of hemodynamic variables. In particular, 
the blood pressure in tumor vessels ranges from typical arterial values to venous values as these 
vessels connect both sides via shunts by a continuously varying pressure potential. Consequently 
the pressure in arteries, over which the tumor as grown, is lower than normal while the pressure 
in former veins is elevated. 

In the flow computations we flx the pressure difference between supply-and drain, which means 
that the flow rates in the major vessels are determined by the total resistance of the network. 
Lowering this resistance by arterio- venous shunts increases the blood throughput, which is evident 
in Fig. 10c when observing the increased flow rates in the thickest vessels. This increase is in 
particular pronounced on the venous side. The flow rate in the initial network grow slightly 
stronger with the radius on the arterial side than on the venous. Therefore the difference between 
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tumor network and normal network is less pronounced for arteries. Interestingly many of the 
dilated capillaries at r'™"^-' fall below the normal average flow rate. 

In contrast, most non-dilated small capillaries which are obviously located in the outer rim fall 
in the same ranges than normal vessels. A few vessels though show significantly lowered flow and 
velocity values. In section [3.41 we study the impact of them upon drug delivery. 

3.3 Parameter dependencies 

Here we discuss the role of the most sensitive parameter dependencies and model changes. The 
large shear force variations among capillaries in the original vasculature leads us to the question 
whether an absolute threshold for the collapse criterion is appropriate. Therefore we also checked 
the result of setting p'-"^ oc 1.0 - {f / finit) / f''"^ for f / finit < f^''\ where finit is the shear force in 
the of the original vasculature. The result is less MVD fluctuations in the samples and also locally 
a more homogeneous distribution. 

In the present model variant, vessels can collapse any time after degradation [w < w'-'^^), 
depending on shear-stress. It is however also an option to restrict vessel collapses close to the 
tumor boundary. For example, it has been suggested that the death of tumor cells releases solid 
pressure from nearby vessels which is otherwise exerted by the tumor (Griffon-Etienne et al., 
1999). In previous papers we assumed the existence of a stable radius, preventing collapses in the 
center where vessel radi are larger. Indeed, the EphB4 signaling mechanism which is related to 
circumferential growth, can also lead to reduced leakiness, tightened EC junctions and increased 
endothelium/pericyte interaction (Erber et al., 2006), making an actual stability improvement 
plausible. Again, we have implemented this very simplistically by modulating the collapse proba- 
bility: p^^' = p^ing = P^fuu ■&{x — xt — S^'^^), where Q is the Heaviside function, xt the radial tumor 
extend, and (J^^-* is the ring width measured from the invasive edge. The result obviously depends 
strongly on the time vessels spend in the collapse region, given implicitly by S^'^\ Aw, the radius 
r and the tumor expansion rate. If this time is of the order of the mean survival time At/p^'^'^ a 
transition occurs toward a highly vascularized center, reminiscent of the "percolation transition" 
analyzed previously (Bartha and Rieger, 2006; Welter et al. 2007; Lee et al., 2007). In this case the 
random collapse-process also plays a more dominant role because many weakly-perfused vessels 

(c) 

can survive whereas, using p^^n , these vessels would eventually regress if one just waits sufficiently 
long. Fig. 11 shows respective simulation results, using = 400Aim, /(c^™'^^) = 3Pa and afl other 
parameters unchanged. One can see increased formation of high-MVD hot spots and less isolated 
threads than for example in Fig. 5. The MVD is slightly higher, but easily tunable by ^(^."laa;) 
j{c,max) _ Generally, p'^^l-^g also leads to less fluctuations with respect to initial networks apparently 
because initially thick vessels provide a backbone which never collapses due io w > w^'^^ while 
being in the "collapse ring" . 

By default sprouting is allowed from veins and capillaries but not from arteries. Enabling 
sprouting from all vessels naturally leads to ca. 40 % higher MVD in the boundary, while the 

(c) 

central MVD also depends on the collapse model. Using parameters with p^ing that the effect 

(c) 

of Prirag becomes significant, i.e. near a transition to fully vascularized tumor, the MVD will 
increase, but otherwise not. Because the boundary becomes more homogeneously vascularized 
also the O2 field becomes more homogeneous. But with the current O2 proliferation threshold 
^(p»"o') 0.9(co) this does not alter the tumors expansion rate due to TCs proliferating around 
smaller I0W-O2 regions and enclosing them. We expect this to be different in reality. TCs would 
likely migrate into I0W-O2 regions driven by the pressure of the surrounding tissue where they 
remain in a quiescent state. 

3.4 Drug flow 

In order to assess the effects of typical tumor network morphology on transport of drugs into and 
through the tumor, we analyze the time-dependent concentration distribution of drug during a 
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continuous injection into the blood stream, starting at arterial boundary nodes. 

The computational model for drug flow in a vascular network is described in detail in (Welter 
et al., 2007; McDougall et al., 2002). The starting point is a given configuration for the vascula- 
ture in our model with precomputed/prescribed variables for flow, flow velocity, segment length, 
and radius of the pipes in the network. In addition, a mass parameter m is now associated with 
segments describing the amount of drug in the segments volume. The mass content m is determin- 
istically updated in successive time steps as follows: First the drug amount flowing out of segments 
is determined and added to respective node-mass variables. Under the assumption of perfect mix- 
ing, the nodal masses are then redistributed into further downstream segments. Thereby mass 
conservation is strictly enforced. A detailed description can be found in (Welter et al., 2007). The 
most severe limitation of this model is that there is no exchange with extra vascular space and 
therefore also no uptake by the tumor. 

The results in the following were obtained with an continuous injection into the vasculatures 
from the case = l. From t=0 on, blood with solute drug at cone. C'(™**) = l flows into 

the vasculature which is initially filled with "clean" blood. We don't consider bolus injections 
because the flow rates arc of the order of mm/s, which is sufficient to saturate 80-100% of the 
vessels within seconds, depending on the network configuration. 

Fig. 15 shows a sequence of snapshots from configuration (a) t=0 .. 6.6s. Drug enters the 
system via the arteries and flows downstream with a sharp transition at the drug/clean interface. 
When vessels merge in upstream direction mixing with clean blood occurs and so the concentration 
increases moderately in the veins. After 20 s near full saturation is achieved. Dilation of tumor 
internal vessels and apparent direct connections to feeding arterioles lead to comparably fast filling 
with drug, whereas a few regions in the highly vascularized boundary take significantly longer to be 
filled. Depending on the exact network configuration we have consistently encountered such strong 
variations in drug delivery to the outer rim. It is likely to be caused by the shunts which lead to 
decreased flow rates in the surrounding tissue. Fig. 10c also shows the existence of capillary sized 
tumor vessels with very low flow rates of the order of lO^/xm ^/s, corresponding to an average flow 
velocity of the order of 10/um/s. The average flow velocity over all vessels in the original network 
is ~ 1.8mm/s. 

Quantitative analysis of drug efficiency impact is only possible to a limited degree due to the 

lack of tumor uptake modeling and because even few surviving TCs can grow a new tumor mass. 
None the less we measured statistics for the total amount of drug in the vasculature. Fig. 16a for 
example, shows the fractional length of the tumor vessel-network for which the maximum drug- 
concentration was larger than indicated on the "c"-axis. After 20 s, on average 94 ±3% of the 
vasculature had been exposed to a drug concentration of at least 80%. For higher concentrations 
the curve decreases drastically. It also shows evidently by the drop at c=0 that 4 ± 3% had not 
been drug perfused at all. In absolute numbers the latter corresponds to 15 ± 12mm. Fig. 16b 
shows the fractional length of the tumor vessel -network for which the concentration was larger 
than the threshold c = 0.25 for a total duration greater than indicated on the x-axis as tg- This 
measurement was also done at t=20s, thus the exposure time cannot be longer than that. One can 
see, for example, that 93 ±4% of the vasculature had been exposed 50% drug concentration longer 
than 10s. Considering the flow rates/velocities mentioned above one can expect the remaining few 
millimeters to be drug perfused within a time frame of minutes. We therefore still believe that the 
geometry and hydrodynamic flow characteristics do not pose an inherent problem to drug delivery. 

4 Hot-spots and spatial inhomogeneities 

Already by visual inspection one observes a coincidence of the location of hot spots in the tumor 
vasculature with the location of the thicker vessels. For example in configuration (a) where an 
artery is close to a vein, or vice versa, an increased number of vessels survive in between. On 
the other configuration similar behavior can be observed on smaller scale. The blood pressure in 
thick vessels is approximately constant compared to the steep drops in the capillary bed. This is 
due to the orders of magnitudes larger flow conductivity. It is plausible that in (a) for example, a 
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high global vertical pressure drop leads to high shear forces in between the two major vessels and 
therefore increased survivability. Furthermore the asymmetry in the shear stress distribution over 
the vascular trees (meaning that given vessels of a certain radius, the shear stress in arteries is 
usually higher than in veins; see Fig. 10), apparently leads to hot spot formation preferably close 
to the arterial branches. 

This motivates a quantitative analysis of the relation between inter-vascular "blood-pressure" 
gradients in the original vasculature and the local MVD of the tumor vasculature. This is done as 
follows: Since new vessels formed during the process of sprouting angiogenesis will experience a 
shear force proportional to the pressure drop from one end to the other. Therefore we first compute 
a pressure field P{x) as solution to the Laplace equation AP = subject to the boundary condition 
that P is set to the blood pressure at vessel sites. The resulting field interpolates pressures between 
adjacent vessels, and would change perfectly linear between two infinite parallel vessels. At any 
given point, we can obtain the magnitude of the gradient , which we use to compute a spatial 
average of ||VP||. We do this for t =0h. The resulting mean value (||V-P||), is plotted together 
with an estimate for the tumor MVD at t =1200h (defined as the fraction of occupied sites) as 
displayed in Fig. 12. The main plot shows the MVD versus (||VP|[) over entire tumors with one 
data point per simulation run. The correlation coefficient is usually very large ~ 0.9 for this global 
measurement. Analogously we analyzed single tumor instances where the averaged is taken locally 
over randomly distributed small regions (discs with 150/im radius). The inlet shows the resulting 
distribution of one such measurement. The correlation coefficient in this instance is already low 
(« 0.4). On average we observed it to vary with the model details and even with the initial 
vasculature between 0.2 and 0.4. 

Fig. 13 is a visualization of the statistical correlation observed in the above analysis. We take 
the data fro the system configuration in the simulation run shown in Fig.5. The bottom image (c) 
shows the site occupation by the vessel network. It is easy to identify the tumor network by its 
typical structure. In (a), where ||VP|| is displayed, one observes that indeed the local gradients 
are stronger in the area close to the major artery in the top half, than close to the vein in the 
lower half. The snapshots on the right column (b) and (d) show respective mean value samples 
of the left hand side. These distributions show the behavior that zones with elevated (||VP||) are 
also probable locations of high MVD. 

A plausible argument for the existence of the correlation between hydrodynamic characteristics 
of the initial network and local morphological features in the evolved tumor vasculature is the 
following: In regions with originally strong inter-vascular pressure drops, tumor vessels would be 
less prone to collapse because sprouting forms new connections whose shear-stress is proportional 
to ||VP||. Thus, shear-stress stabilization would prevent regression in high-||VP|| regions. Since 
single collapse events lead to long-ranged collapses of adjacent network sections, we believe that 
therefore local measurements hardly show correlations. Note that the intra- vascular blood-pressure 
at in-and outflow vessels is prescribed and depends on the vessel radius. Therefore only the 
positioning of initial vessels leads to the variety of emerging tumor morphologies. 

We further quantify the spatial inhomogeneities by probability distributions for local MVD 
(Fig. 14a), necrotic region size (b) and hot spot size (c). The distributions are estimated via 
histograms of the respective local measures, whereby we insert data from multiple (40) runs 
into the same histogram in order to get a smooth curve with reasonably small bin. The MVD 
in Fig. 14a is estimated as fraction of occupied sites within boxes (250/Ltm in size) of a regular 
grid. The plot displays two peaks. One large at low MVD«0.07 which decays algebraically with 
exponent ^ 1.4 till smaller peak at MVDwO.45. The latter stems from the high-MVD zone in the 
tumor periphery. Not included is a peak which appears at MVD=0 due to large non-vascularized 
regions. The necrotic-region size distribution in Fig. 14b plots the probability to find a connected 
cluster of dead TCs with a certain volume, i.e. number of cells. The hot-spot size distribution 
in Fig. 14c plots the probability to find a connected cluster of a certain number of sites where the 
local MVD is larger than a threshold value (here 0.15). Therefor we compute an estimate for the 
MVD for each lattice site as the fraction of occupied sites within a 250/im diameter disc. Both of 
the latter distributions exhibit purely algebraic decay also with the exponent ~ 1.4. 

The fractal dimension d/ is often used to quantify tumor networks and has been found capable 
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to distinguish them from normal vascular networks (Baish and Jain, 2000). Therefore we do box 
counting analysis, which is usual practice to estimate df for natural objects. Perfect fractals show 
a power law N oc e"**^ in the number of touched boxes A'' vs. size e, yielding a constant slope in a 
log-log plot. Wo typically observe a parabolic shape of the local slope curve with a small constant 
plateau over half a decade. Due to the problems associated with df estimates, such plateaus are 
often considered sufficient to extract df, see the discussion in (Chung and Chung, 2001). In this 
context one also speaks of box counting dimension. On average wo obtain df = 1.76 ± 0.03 for the 
final internal tumor network; df = 1.79 ± 0.03 including high-MVD peripheral region. The error 
is given by the root mean square deviations among 40 simulation runs with the same parameter 
set. Fluctuations of this magnitude are very significant, which is why we consider these results as 
unreliable and also show no data here. Furthermore resembling previous results, we still find our 
df estimate being primarily determined by the MVD. 

5 Discussion 

The conclusions that can be drawn from the theoretical model for vascular remodeling of an 
arterio-venous network during tumor growth are manifold: As has already been conjectured on 
the basis of an analogous model with a regular or lattice-like blood vessel network (Welter et al., 
2006) the global characteristics of the emerging tumor vasculature is expected to be independent of 
the specific details of the initial vasculature: the resulting morphology is compartmentalized into 
a highly vascularized tumor perimeter, a tumor periphery with large vessels density and dilated 
vessels and a central region containing necrotic regions with a low micro- vascular density threaded 
by extremely dilated vessels. The basic mechanisms leading to this morphology are identified 
as sprouting angiogenesis in the tumor perimeter, a switch of the vascularization program to 
circumferential growth plus vessel regression within the tumor. 

On top of these global feature our model predicts spatial inhomogeneities, visible in drastic 
MVD variations. These are strongly correlated with the local characteristics of the initial arterio- 
venous vessel network, which are absent in a regular homogeneous initial vasculature. In simulation 
runs with identical parameters but different arterio-venous network configuration, cases can be 
observed ranging from tumors threaded mostly by individual vessels, to tumors that exhibit many 
dense clusters connected by few short chains. The reason for this variability the asymmetry of the 
shear force distribution between the arterial and venous side of the vascular trees. 

Depending on the details of the initial network in our model the tumor vasculature develops 
isolated highly vascularized clusters connected by thick vessels. These "hot spots" are also observed 
in real tumor and serve as an important diagnostic tool in cancer therapy. Already by a visual 
comparison of the starting network with the final tumor network one observes that hot spots form 
more frequently in those regions where the starting network contained predominantly arteries. We 
therefore analyzed the correlation between various local hydrodynamic quantities of the original 
network (blood pressure, blood pressure gradient, blood fiow) and the local MVD in the tumor 
network and found a significant correlation between local blood pressure gradient of the original 
arterio-venous network and the most probable locations of hot spots in the tumor vasculature. This 
is plausible since a high pressure gradient within the vessels implies a high shear force exerted by 
the blood on the vessel walls, which stabilizes the vessel against collapse and therefore leads to an 
increased vessel survival probability when the tumor has grown over this region. Since we expect 
these mechanisms also to be at work in real tumors it should therefore in principle be possible to 
predict on the basis of an analysis of the blood vessel network in the healthy tissue, where most 
likely the tumor vasculature develops its hot spots. We note that in our model this hydrodynamic 
mechanism plays a dominant role in the the hot spot formation, whereas potential local variations 
in pro- and anti-angiogenic effectors within the tumor were not involved but could also play a role 
in this process. 

Remarkably an association between the MVD of the original vasculature in normal tissue and 
hot spot vascular density in the tumor has been reported for carcinoma (Hockel et al., 2001). 
Here the MVD of normal tissue was measured far away from the tumor and it was found that it 
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correlates well with the MVD of the hot spots in the tumor. Our interpretation of this result is 
that sprouting angiogenesis is not strong in these tumor samples (in fact the MVD in the tumor 
periphery was increased by only ca. 20%) and hot spots in the tumor represent regions of the 
original vasculature that could survive vessel collapse due to locally increased pressure gradients. 
For a reliable check of our hypothesis one would have to analyze the pressure gradients in the 
vasculature of the normal tissue before the tumor grows over it, i.e. under experimental conditions 
of an implanted tumor for instance. 

Various probability distributions that quantify the spatial inhomogeneities of the tumor net- 
work in our model turn out to display an interesting behavior, too: The probability distribution 
of local MVD values as well as the probability distribution of the hot spot volume (defined as the 
size of connected clusters of regions with MVD larger than the original MVD) have an algebraic 
tail with exponent ~ 1.4. Also the size distribution of the connected clusters of necrotic tumor 
regions showed an algebraic tail (again with exponent ^ 1.4). These power laws indicate a critical 
state in accordance with the fractal properties of the tumor network: It implies the absence of a 
particular lengths scales over which size distribution would decay exponentially and resemble the 
size distribution that one encounters in percolation at the critical point, the percolation threshold. 

When analyzing the blood flow characteristics in emerging tumor vasculature of our model 
we found an interesting phenomenon which is naturally absent in models with a regular starting 
configuration: Thick arterioles and vcnoules provide a well conducting support structure around 
the tumor. Since the total pressure difference between the tree roots is fixed, the transported 
blood volume is given by the total flow resistance of entire vascular tree. Dilation of a few selected 
vessels forming a path between the tree roots can remove a bottlenecks formed by thinner vessels. 
Angiogenesis which provides additional vessels thereby facilitates the creation arterio- venous short- 
cuts, or shunts, through multiple partly disjoint paths. This leads to an increased blood flow 
through the tumor vasculature when compared with the starting vasculature. In contrast to this 
in a regular network the total flow resistance is dominated by the network outside of the tumor 
(Bartha and Rieger, 2006, Welter et al. 2007) and the flow cannot not increase via the dilation 
of tumor internal vessels. There arc experimental studies which agree with the predictions of our 
model. In (Sahani et al., 2005) perfusion parameters in rectal cancer were measured via computer 
tomography where consistently increased blood flow is reported by approximately a factor of two 
compared to normal tissue. There it is also argue that angiogenesis facilitates the creation of 
arterio- venous shunts which bypass the capillary network - the same mechanism by which blood 
increases in our model. 

These artcrio-venous shunts should also lead to an increased blood pressure range in tumor 
vasculatures: They connect vessels with high pressure to vessels with low pressure, which then 
varies continuously within the connecting vessel. Indeed we flnd a wide spread among pressures 
over a short range of radii near r^ax in our model, which is somewhat in contrast with the 
prevailing view that the pressure as function of vessel radius is reduced in the venous part in 
tumors (Jain, 1988). Our result indicate that the measure used to assess the relation of tumor 
blood flow to normal tissue has critical impact on the results. It has been shown that measuring 
perfusion as function the spatial location can give rise to elevated flow rates, whereas comparison 
of flow rates between vessels of a certain radius may show either increased as well as decreased 
flow rates depending on the radius and the location in the vascular hierarchy. 

Our drug flow computations again showed that a drug bolus, which is injected into the source 
arteries of the blood vessel network in our model for some time, reaches all perfused blood vessels, 
although the thin capillaries in the highly vascularized tumor boundary needs a longer time to 
become filled with drug. Blood-borne delivery of therapeutics into the tumor- vascualture does not 
appear to be an obstacle for a successful chemotherapy. The reasons for failure of drug delivery 
to tumor cells arc more likely related to drug uptake the drug transport through the tumor tissue 
(Minchinton and Tannock, 2006). Since blood is in contact with the interstitial fluid, due to leaks in 
tumor vessels, fluid exchange can take place with the ECM. This has been shown to lead to steady 
states with elevated IFP levels, where the difference between micro-vascular pressure (MPD) and 
IFP generally lower than normal (Boucher et al., 1996). Since convective transport is driven by 
pressure differences, high IFP could pose a barrier to drug delivery (Hassid et al., 2006)). On 
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the other hand, leakiness and MVP-IFP gradients could lead to premature release predominantly 
in locally restricted regions around vessels where blood enters the tumor. Vessels in the outflow 
regions would thus be depleted of drug. Locally released drug would then be transported by IFP 
gradients out of the tumor. Furthermore drugs usually consist of large macromoleculcs. Their 
low diffusibility through the vessel wall and generally lower diffusibility than O2 could lead to 
situations where sufficient O2 reaches certain TCs to let them remain viable, but not enough drug 
reaches them to kill them off due to the lower diffusion radius. 

In future work it would be useful to study the remodeling of a three dimensional arterio- 
venous network. Getting access to a realistic initial network is likely to pose a hard problem. In 
particular if it must be artificially c;reated. Otherwise the extension to 3d is straight forward and 
has already been done in (Lee et al., 2006) with similar results. Further aspects to assess are 
interstitial fluid transport and O2 release with associated intra-vascular O2 decrease. Furthermore 
implementation of a tumor model which supports quantitative assessment of stresses in the tissue 
appears mandatory for studies of irregular tumor morphologies. 
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Tables 

Table 1: Listing of model parameters with values for our base case simulation runs. 



Parameter 


Value 


Description 


Reference 




60 /xm 


Lattice const. (Tree-generator) 




5 


10 iJ.m 


Lattice const. (Tumor-model) 




L 


1200 


Lattice size (Tumor-model) 




A/ 


1 h 


Time stc]) 




in* =0)1 


1000 


Initial number of TCs 






0.02 


O2 source coefficient 


(Secombet al. 2004) 




(80 iim)-'^ 


consumption in normal tissue 


(Carmeliet and Jain, 2000) 






consumption tumor tissue 




Q(prol) 


0.3 « 0.9(co) 


TC O2 prol. threshold 






0.45 


dimensionless blood-02 level 




fi{death) 




TC hypoxia threshold 






200 Aim 


Growthfactor diffusion radius 


(Nehls et al. 1998) 




10-4 


sprouting GF threshold 




(uo) 
TC 


100 h 


TC survival time under hypoxia 


(Yu et al. 2002) 


Aprol) 


10 h 


TC proliferation time 




.{sprout) 
EC 


5 h 


Sprout initiation/extension time 




(dii) 

''EC 


40 h 


Vessel dilation time 


(Dome et al. 2002) 




100 h 


Time till sprout regression 


(Nehls et al. 1998) 


.{switch) 
'■EC 


25 h 


Time till switch to circumferential growth 




Aw 


0.04 /im/h 


Wall thickness/stability decrease rate 




j,{init) 


4 yum 


Initial vessel radius 






25 /xm 


Max. vessel radius 


(Dome et al. 2002) 


lispr) 


20 urn 


Min. distance between junctions 


(Dome et al. 2002) 


j{c,max) 


1 Pa 


Peak critical shear force 




p{c,max) 


0.1 


Peak collapse probability 
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Figure captions 



Fig.l: Illustration of the model state: All discrete elements are aligned at the triangular lattice 

shown in the background. Blue bars represent vessels occupying bonds. The set containing all 
vessels is denoted by V. One vessel is highlighted by a black frame, it's junctions to other vessels 
are denoted by xi.X2- Tumor cells (TCs) can occupy lattice sites with a one to one relationship. 
The set containing TCs is denoted as T. One TO is displayed as yellow hexagon. 

Fig. 2: Sketch of the dynamical processes in the model: (a) TC Proliferation, (b) TC Death, 
(c) Sprout formation, (d) Sprout migration, (e) Vessel regression and removal, (f) Dilation and 
wall degradation. See text for details. 

Fig. 3: Examples of starting configurations of the model at t =0h. We chose four different 
basic layouts for the vascular network denoted as A,B,C,D. They differ in the placement and 
lengths of major supply-and drain vessels, which are randomly distributed. This is indicated in 
the top right inlets. Arteries are red, veins blue, black bars show intervals from which starting 
locations are drawn while the wedge shaped part of the red/blue bars indicate intervals from which 
starting lengths are drawn. The color code of the resulting networks shows the blood pressure 
distribution from 2kPa (red) to 12kPa (blue). The small yellow area in the system center is the 
tumor nucleus initialized by eden growth to maximum of 1000 TCs. The background in normal 
tissue/extracellular matrix is colored green if the growthfactor is sufficient to enable sprouting. 
Otherwise those areas are white. 

Fig. 4: Snapshots of the dynamical evolution at times t =0, 200, 600h of the model in a single 
simulation run starting in configuration A (shown in Fig. 3a) at t = 0. Viable TCs are bright 
yellow, while necrotic regions are dark yellow. Vessels are color coded by blood pressure as in 
Fig. 3. The background in normal tissue/extracellular matrix is colored green if the growthfactor 
is sufBcient to enable sprouting. Otherwise those areas are white. 

Fig. 5: Final configuration at t =1200h for the simulation run depicted in Fig. 4. The color 

code is identical to Fig. 4. Note how vessels in the tumor center serve as arteriovcneous shunts. 
They have dilated diameters and thus carry high blood throughput. The formation of hot spots 
with increased local MVD is also evident. 

Fig. 6: Final configurations at {t =1200h) of simulation runs of the model startig with the 
initial configurations B, C, D shown in Fig. 3. The color code is identical to Fig. 4. The formation 

of hot spots is clearly visible in the central region in B, but not so in C and D. There increased 
MVD can be observed in a few locations close to the tumor boundary. 

Fig. 7: Afagnified view of the tumor boundary showing in detail how the vasculature is remod- 
eled close to the invasive edge. Alternating with very high MVD on small scale, there are "holes" 
of the size of the intervascular distance in normal tissue. This is an effect of prohibiting sprouts 
from arteries, which would otherwise fill the holes, leading to a homogeneous vascular density. 
Behind the peripheral zone, the MVD drops rapidly. One can further see a few isolated dilated 
vessels directly connected to arterioles/ venoules in healthy tissue beyond the invasive edge. For 
the color code see the caption of Fig. 4. 

Fig. 8: Flow rate distribution in the tumor vessel network depicted in Fig. 5. The color code is 
on a logarithmic scale as indicated on the left side. The flow rate here is the blood volume per time 
through the vessel cross-section, closely related to the blood flow per tissue volume as commonly 
measured in experiments. Dilated tumor-internal vessels usually exhibit high flow rates. If one 
would measure in/outflow per tissue volume element as in dynamic MRT (magnetic resonance 
tomography) measurements (Pahernik et al., 2001) one would observe hot-spots located at the 
respective high-MVD zones. 
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Fig. 9: Radial distributions of various quantities characterizing the dynamical development of 
the compartmentalization of the tumor vasculature in the model: The radius r here denotes the 
distance from the tumor center. One data point represents an average over a concentric 50/Ltm wide 
shell, and over 40 simulation runs, 10 per configuration A,B,C,D (see text or Fig. 3). (a) shows 
the micro vascular density, (b) the O2 level and (c) the tumor cell density. The right column 
shows blood flow related variables: (d) the vessel radius, (e) the flow rate and (f) the shear force, 
which are in contrast to (a),(b),(c) given as the average over actually occupied sites within a shell, 
instead of over all sites. For t =1200h also the local variations (root mean square deviations) are 
indicated as error bars. 

Fig. 10: Scatter plots showing hydrodynamic quantities for individual vessel segments as a 
function of the vessel radius: (a) blood pressure, (b) wall shear force, and (c) flow rate through 
the vessel cross section. Each data point is a sample from a randomly chosen position, uniformly 
distributed over the vascular network. Samples from both, the initial vasculature at t =0h, and 
the tumor tumor network at t =1000h are shown. For the former, arteries are displayed in red 
with negative radius, veins arc blue and capillaries are pink. Tumor vessels that were initially 
arteries are displayed yellow, those that were veins are light blue. New vessels and capillaries are 
black. Dilation causes the latter vessels to span radius values from 5 to 25/xm (the maximum 
dilation radius r*^™"^)), but due to the lack of hierarchical organization in the tumor they cannot 
be contribute to arteries or veins. Therefore respective data points were put with probability 1/2 
on one or the other side. Note that tumor vessels capped at r^™"^^ display a much larger range of 
values than normal (see text). Also note that (b) and (c) are presented as linear- log plots. 

Fig. 11: Final configuration at t =1200h of a simulation run using an alternate vessel collapse 
rule where vessel collapses are restricted to a thin band (here 400/um in width) behind the invasive 

(c) 

edge. The collapse probability p^ing is thus modulated with an appropriate term which depends 
on the distance to the tumor center (see text). In this instance, parameters are such that the mean 
survival time of instable vessels is of the order of the time they spend in the "collapse region" , 
leading to a dominant role of the random collapse process for the resulting morphologies, and 
putting the system near a transition to a fully vascularized tumor. The parameter set is the same 
as in the base case, except for the critical shear force /(c-™"^^) which is set to 3Pa, instead of IPa. 

Fig. 12: Local MVD in the tumor vasculature versus local pressure gradient (|AP|) (i.e. blood 
pressure differences between neighboring vessels) in the starting vasculature. In the main plot 
data points are averaged over different runs of the base case scanrio, the averages are taken over 
the entire tumor interiors (excluding the vascularized boundary region) at t =1200h. The inlet is 
a scatter plot for a single simualtion run, where the an average is performed locally over a small 
disc with radius 150/xm (see Fig. 13). 

Fig. 13: Digitized data processing of network configurations demonstrating the correlation 
between the intervascular pressure gradient in the starting network with the local MVD in the 
tumor network, (a) Local magnitude of the "intervascular pressure" gradient || AP|| in the initial 
network (see text), (b) Smoothing of the data in a: locally averaged values of || AP|| were computed 
at randomly distributed locations, each by averaging over a 150/zm-radius disc, then the space 
closest to some data point is filled with a gray scale value proportional to (|| AP||). (c) Local MVD 
of the tumor network, (d) Smoothing of the data in c, analogous to b, yielding an estimate for 
locally averaged MVD. Value ranges are 0-45 (black to white) for ||AP|| and 0-0.4 for the MVD. 
Note that in d the local MVD is highest (brightest) in regions in the upper half of the network, 
and in b the local precussure gradient is highest (brightest) also in regions in the upper half of the 
network. 

Fig. 14: Shows probability distributions in log-log plots, for (a) the local MVD, given as the 
local average over 250/zm wide boxes, (b) the volume of necrotic tissue clusters, defined as the 
number of sites in respective connected components of dead tissue, (c) the volume of vessel 
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hot-spot areas, defined as the connected components of regions where the local MVD exceeds a 
threshold (0.15). The distributions are generated by binning observed values in histograms over 40 
simulation runs (at t =1200h). Note that the distributions show algebraic decay. In this instance 
in particular with the same exponent within the error bounds which are of the order of 2%. 

Fig. 15: Snapshots of the drug-flow simulation using the configuration A tumor at t =1000h. 

The color code shows the drug concentration as indicated on the scale bar. Note that drug reaches 
most parts of the vasculature quickly in a few seconds at the maximum concentration. In the 
T=6.6s snapshots a few vessels in the boundary region are not yet perfused. The time until 
drug saturates the complete vasculature is of the order of minutes. We can observe this behavior 
generally for all resulting networks. 

Fig. 16: (a) shows the amount of vessels given as the fraction of total tumor network length 
that has been exposed to a drug concentration larger than the indicated concentration Cmax (at 
t=20s of the drug simulation). The black line represents the average over 40 runs with different 
networks, (b) shows the amount of vessels givc^n as the fraction of total tumor network length 
that has been exposed to a drug concentration larger than 0.25 for longer than indicated time tg- 
The black line represents the average over 40 runs with different networks. Like in Fig. 15, tumor 
networks from the base case at t =1000h where used. 
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